# after install the packages below, run and load them
  library(readxl)

  # install.packages("remotes")
  # remotes::install_github("mdelacre/deffectsize")
  library(deffectsize)
  library(statsExpressions)
  
  library(effectsize)
  # devtools::install_github("tomfaulkenberry/anovaBFcalc")
  # anovaBFcalc::anovaBFcalc()
  library(anovaBFcalc)
  library(BayesFactor) 
  
  library(ARTool)
  library(dplyr)
  library(rlang)
  
# update R
  # library(installr)
  # updateR() 
  # update.packages(checkBuilt=TRUE, ask=FALSE)
  # installed.packages()[,c('Package','Version','LibPath')]
--------------------------------------------------------------------------------
                       Participant characteristics
--------------------------------------------------------------------------------
# Effect Size (hedges g*) using package "deffectsize"(="effectsize"or"statsExpressions")
  # Age 
  cohen_CI(m1=22.0833,m2=22.3056,sd1=2.1696,sd2=1.9540,n1=36,n2=36,
           conf.level=0.95,var.equal=FALSE,unbiased=TRUE,alternative="two.sided")
  # hedges_g(Age~Group, data=participant, paired= FALSE, pooled_sd=FALSE) # using "effectsize"
  # two_sample_test(x=Group, y=Age, data=participant, type="parametric", var.equal=FALSE, effsize.type="g") using "statsExpressions"
  
  # Education 
  cohen_CI(m1=15.9444,m2=16.1111,sd1=2.3415,sd2=2.0358,n1=36,n2=36,
           conf.level=0.95,var.equal=FALSE,unbiased=TRUE,alternative="two.sided")
  
  # MaRs-IB
  cohen_CI(m1=77.5310,m2=77.4171,sd1=10.9490,sd2=11.4466,n1=36,n2=36,
           conf.level=0.95,var.equal=FALSE,unbiased=TRUE,alternative="two.sided")
  
  # Handedness
  cohen_CI(m1=81.0222,m2=83.2639,sd1=16.8924,sd2=15.1902,n1=36,n2=36,
           conf.level=0.95,var.equal=FALSE,unbiased=TRUE,alternative="two.sided")
  
  # Years of music training using package "statsExpressions"
  participant <- read_excel("C:/Users/ALIENWARE/Desktop/participant.xlsx")
  names <- c('Subject','Group')
  participant[,names] <- lapply(participant[,names],factor)
  str(participant)
  summary(participant)
  
  one_sample_test(data=participant,x=Music_Years,test.value=0,
                  type="parametric",k=3L,effsize.type="g")
  # hedges_g(Music_Years~ 1, data=participant, mu=0) # using "effectsize"
 
  # Musical Training
  cohen_CI(m1=5.56349,m2=1.31349,sd1=0.46567,sd2=0.25053,n1=36,n2=36,
           conf.level=0.95,var.equal=FALSE,unbiased=TRUE,alternative="two.sided")
  
  # Perceptual Abilities
  cohen_CI(m1=5.98148,m2=3.72531,sd1=0.67010,sd2=0.94896,n1=36,n2=36,
           conf.level=0.95,var.equal=FALSE,unbiased=TRUE,alternative="two.sided")
  
  # Singing Abilities
  cohen_CI(m1=5.32937,m2=2.66667,sd1=0.74099,sd2=1.06795,n1=36,n2=36,
           conf.level=0.95,var.equal=FALSE,unbiased=TRUE,alternative="two.sided")
  
  # Active Engagement
  cohen_CI(m1=5.24383,m2=3.09568,sd1=0.67707,sd2=1.01840,n1=36,n2=36,
           conf.level=0.95,var.equal=FALSE,unbiased=TRUE,alternative="two.sided")

  # Emotions
  cohen_CI(m1=5.75926,m2=4.34722,sd1=0.72860,sd2=0.97376,n1=36,n2=36,
           conf.level=0.95,var.equal=FALSE,unbiased=TRUE,alternative="two.sided")
  
  # General Sophistication
  cohen_CI(m1=5.39815,m2=2.58179,sd1=0.50614,sd2=0.76756,n1=36,n2=36,
           conf.level=0.95,var.equal=FALSE,unbiased=TRUE,alternative="two.sided")
  
--------------------------------------------------------------------------------
                                Stimuli 
--------------------------------------------------------------------------------  
## Recognition rates
  # Effect Size for one-way ANOVA using package "effectsize"
  print(F_to_omega2(5.87, 2.83, 56.70,ci=0.95,alternative="two.sided"),digits=3) # Emotion
  
  interpret_omega_squared(.185, rules="field2013") #This rules are same as ones suggested by Kirk (1996).   
  
  # Bayes Fators (BF) for main effect using "anovaBFcalc" MWS method with alpha = -0.5 (default)
  bf_mws(F=5.87, df1=2.83, df2=56.70, report.as="BF10",alpha=-0.5)  # Emotion 
  
  # Effect Size for Post comparison between emotion
  t_to_d(t= -4.28,df_error=20,paired=TRUE,ci=0.95,alternative="two.sided") # Happiness-Sadness 
  t_to_d(t= -1.54,df_error=20,paired=TRUE,ci=0.95,alternative="two.sided") # Happiness-Fear
  t_to_d(t= -2.97,df_error=20,paired=TRUE,ci=0.95,alternative="two.sided") # Happiness-Anger
  t_to_d(t= -4.57,df_error=20,paired=TRUE,ci=0.95,alternative="two.sided") # Happiness-Neutrality  
  t_to_d(t=  1.75,df_error=20,paired=TRUE,ci=0.95,alternative="two.sided") # Sadness-Fear
  t_to_d(t=  0.98,df_error=20,paired=TRUE,ci=0.95,alternative="two.sided") # Sadness-Anger
  t_to_d(t= -0.36,df_error=20,paired=TRUE,ci=0.95,alternative="two.sided") # Sadness-Neutrality
  t_to_d(t= -1.23,df_error=20,paired=TRUE,ci=0.95,alternative="two.sided") # Fear-Anger
  t_to_d(t= -2.48,df_error=20,paired=TRUE,ci=0.95,alternative="two.sided") # Fear-Neutrality
  t_to_d(t= -0.75,df_error=20,paired=TRUE,ci=0.95,alternative="two.sided") # Anger-Neutrality
  
  
  interpret_cohens_d(c(−0.96,-0.34,-0.66,-1.02,0.39,0.22,-0.08,-0.28,-0.55,-0.17),rules="cohen1988")
  
  # BF for paired t-test using "BayesFactor" 
  exp(ttest.tstat(t= -4.28, n1=21,rscale=0.707)[['bf']]) # Happiness-Sadness  
  exp(ttest.tstat(t= -1.54, n1=21,rscale=0.707)[['bf']]) # Happiness-Fear 
  exp(ttest.tstat(t= -2.97, n1=21,rscale=0.707)[['bf']]) # Happiness-Anger
  exp(ttest.tstat(t= -4.57, n1=21,rscale=0.707)[['bf']]) # Happiness-Neutrality  
  exp(ttest.tstat(t=  1.75, n1=21,rscale=0.707)[['bf']]) # Sadness-Fear   
  exp(ttest.tstat(t=  0.98, n1=21,rscale=0.707)[['bf']]) # Sadness-Anger  
  exp(ttest.tstat(t= -0.36, n1=21,rscale=0.707)[['bf']]) # Sadness-Neutrality
  exp(ttest.tstat(t= -1.23, n1=21,rscale=0.707)[['bf']]) # Fear-Anger  
  exp(ttest.tstat(t= -2.48, n1=21,rscale=0.707)[['bf']]) # Fear-Neutrality 
  exp(ttest.tstat(t= -0.75, n1=21,rscale=0.707)[['bf']]) # Anger-NNeutrality 
  
## Intensity ratings
  # Effect Size for one-way ANOVA using package "effectsize"
  print(F_to_omega2(14.50, 2.11, 42.20,ci=0.95,alternative="two.sided"),digits=3) # Emotion
 
  interpret_omega_squared(.386, rules="field2013")  
  
  # BF for main effect using "anovaBFcalc" MWS method with alpha = -0.5 (default)
  bf_mws(F=14.50, df1=2.11, df2=42.20, report.as="BF10",alpha=-0.5)  # Emotion
  
  # Effect Size for Post comparison between emotion 
  t_to_d(t= -4.99,df_error=20,paired=TRUE,ci=0.95,alternative="two.sided") # Happiness-Sadness 
  t_to_d(t= -2.99,df_error=20,paired=TRUE,ci=0.95,alternative="two.sided") # Happiness-Fear
  t_to_d(t= -4.66,df_error=20,paired=TRUE,ci=0.95,alternative="two.sided") # Happiness-Anger
  t_to_d(t=  3.36,df_error=20,paired=TRUE,ci=0.95,alternative="two.sided") # Sadness-Fear
  t_to_d(t=  2.30,df_error=20,paired=TRUE,ci=0.95,alternative="two.sided") # Sadness-Anger
  t_to_d(t= -1.14,df_error=20,paired=TRUE,ci=0.95,alternative="two.sided") # Fear-Anger 
  
  interpret_cohens_d(c(-1.12,-0.67,-1.04,0.75,0.51,-0.25),rules="cohen1988")
  
  # BF for paired t-test using "BayesFactor" 
  exp(ttest.tstat(t= -4.99, n1=21,rscale=0.707)[['bf']])  # Happiness-Sadness  
  exp(ttest.tstat(t= -2.99, n1=21,rscale=0.707)[['bf']])  # Happiness-Fear 
  exp(ttest.tstat(t= -4.66, n1=21,rscale=0.707)[['bf']])  # Happiness-Anger
  exp(ttest.tstat(t=  3.36, n1=21,rscale=0.707)[['bf']])  # Sadness-Fear   
  exp(ttest.tstat(t=  2.30, n1=21,rscale=0.707)[['bf']])  # Sadness-Anger  
  exp(ttest.tstat(t= -1.14, n1=21,rscale=0.707)[['bf']])  # Fear-Anger  
   
--------------------------------------------------------------------------------
                                  Results
--------------------------------------------------------------------------------
## Nonparametric Type III two-way mixed ANOVA on Accuracy using package "ARTool"
  behavior <- read_excel("C:/Users/ALIENWARE/Desktop/data-2SD.xlsx")
  names <- c('Subject','Group','Emotion')
  behavior[,names] <- lapply(behavior[,names] , factor)
  str(behavior)
  summary(behavior)
  
  m <- art(Accuracy~Emotion*Group+(1|Subject), data=behavior)
  summary(m)
  anova(m, response="aligned")
  b <- anova(m, response="art",type=3)
  print(b, verbose=TRUE,digits=3) 
  
  # Effect Size for two-way mixed ANOVA using package "effectsize"
  print(F_to_omega2(62.10,1, 70,ci=0.95,alternative="two.sided"),digits=3) # Group
  print(F_to_omega2(29.04,4,280,ci=0.95,alternative="two.sided"),digits=3) # Emotion
  print(F_to_omega2( 5.67,4,280,ci=0.95,alternative="two.sided"),digits=3) # Emotion*Group
  
  interpret_omega_squared(c(.459, .282, .062), rules="field2013") 
  
  # BF for main effects and interactions using "anovaBFcalc" MWS method with alpha = -0.5 (default)
  bf_mws(F= 62.10,df1=1,df2= 70,report.as="BF10",alpha=-0.5) # Group
  bf_mws(F= 29.04,df1=4,df2=280,report.as="BF10",alpha=-0.5) # Emotion 
  bf_mws(F=  5.67,df1=4,df2=280,report.as="BF10",alpha=-0.5) # Emotion*Group 
    
  ## Nonparametric Type III two-way mixed ANOVA on Response Time using package "ARTool"
  behavior <- read_excel("C:/Users/ALIENWARE/Desktop/data-2SD.xlsx")
  names <- c('Subject','Group','Emotion')
  behavior[,names] <- lapply(behavior[,names] , factor)
  str(behavior)
  summary(behavior)
  
  m <- art(Response_Time~Emotion*Group+(1|Subject), data=behavior)
  summary(m)
  anova(m, response="aligned")
  b <- anova(m, response="art",type=3)
  print(b, verbose=TRUE,digits=3) 
  
  # Effect Size for two-way mixed ANOVA using package "effectsize"
  print(F_to_omega2(14.33,1, 70,ci=0.95,alternative="two.sided"),digits=3) # Group
  print(F_to_omega2(26.59,4,280,ci=0.95,alternative="two.sided"),digits=3) # Emotion
  print(F_to_omega2( 1.51,4,280,ci=0.95,alternative="two.sided"),digits=3) # Emotion*Group
  
  interpret_omega_squared(c(.156, .264, .007), rules="field2013") 
  
  # BF for main effects and interactions using "anovaBFcalc" MWS method with alpha = -0.5 (default)
  bf_mws(F= 14.33,df1=1,df2= 70,report.as="BF10",alpha=-0.5) # Group
  bf_mws(F= 26.59,df1=4,df2=280,report.as="BF10",alpha=-0.5) # Emotion 
  bf_mws(F=  1.51,df1=4,df2=280,report.as="BF10",alpha=-0.5) # Emotion*Group 
    
--------------------------------------------------------------------------------    
## Repeated Measures Correlation
  library(rmcorr) 
  library(ggplot2)
  library(pals)
  library(cowplot)

  # Calculate rmcorr for the combined group 
  behavior <- read_excel("C:/Users/ALIENWARE/Desktop/data-2SD.xlsx")
  names <- c("Subject", "Group", "Emotion")
  behavior[,names] <- lapply(behavior[,names], factor)
  str(behavior)
  summary(behavior)
  
  my.rmc<-rmcorr(participant= Subject, 
                 measure1= Response_Time,
                 measure2= Accuracy,
                 dataset= behavior,
                 CI.level= 0.95,
                 CIs= "bootstrap",
                 nreps= 10000,
                 bstrap.out= FALSE) 
  my.rmc
  
  # Calculate rmcorr for the music group 
  behavior1 <- read_excel("C:/Users/ALIENWARE/Desktop/data-2SD1.xlsx")
  names <- c("Subject", "Group", "Emotion")
  behavior1[,names] <- lapply(behavior1[,names] , factor)
  str(behavior1)
  summary(behavior1)
  
  my.rmc1<-rmcorr(participant= Subject,
                  measure1= Response_Time,
                  measure2= Accuracy,
                  dataset= behavior1,
                  CI.level= 0.95,
                  CIs= "bootstrap",
                  nreps= 10000,
                  bstrap.out= FALSE) 
  my.rmc1
  
  # Calculate rmcorr for the nonmusic group 
  behavior2 <- read_excel("C:/Users/ALIENWARE/Desktop/data-2SD2.xlsx")
  behavior2 <- read_excel("data-2SD2.xlsx")
  names <- c("Subject", "Group", "Emotion")
  behavior2[,names] <- lapply(behavior2[,names] , factor)
  str(behavior2)
  summary(behavior2)
  
  my.rmc2<-rmcorr(participant= Subject,
                  measure1= Response_Time,
                  measure2= Accuracy,
                  dataset= behavior2,
                  CI.level= 0.95,
                  CIs= "bootstrap",
                  nreps= 10000,
                  bstrap.out= FALSE) 
  my.rmc2

## Create plots
  # Get sample size
  n  <- length(unique(behavior$Subject))
  n1 <- length(unique(behavior1$Subject))
  n2 <- length(unique(behavior2$Subject))
  
  tiff(filename="Fig3.tiff", width=20, height=15, units="in", res=800)
  
  p <-ggplot(behavior, aes(x= Response_Time, y= Accuracy, group= Subject, 
                           color= Subject)) +  
    geom_point(aes(colour= Subject)) +
    geom_line(aes(y= my.rmc$model$fitted.values), linetype= 1) +
    ggtitle("Combined Group") +
    ylab("Accuracy (%)") +
    xlab("Response Time (s)") +
    theme_grey() +
    scale_shape_identity() +
    theme(legend.position= "none",
          plot.title= element_text(face="bold", size= 20, hjust= 0.5),
          axis.title= element_text(face="bold", colour="black", size= 16),
          axis.text= element_text(face="bold",size=14),
          panel.grid.minor= element_blank(),
          axis.ticks.length= unit(0.15,"cm")) +
    scale_x_continuous(expand=expansion(add= 0.1),limits=c(0, 4),
                       breaks=seq(0, 4, by=1))+
    scale_y_continuous(expand=expansion(mult= c(0.01,0.01)),limits=c(0, 120),
                       breaks=seq(0, 120, by=30))+
    scale_colour_manual(values= kovesi.rainbow(n))
  
  p1 <-ggplot(behavior1, aes(x= Response_Time, y= Accuracy, group= Subject, 
                             color= Subject)) +  
    geom_point(aes(colour= Subject)) +
    geom_line(aes(y= my.rmc1$model$fitted.values), linetype= 1) +
    ggtitle("Music Group") +
    ylab("Accuracy (%)") +
    xlab("Response Time (s)") +
    theme_grey() +
    scale_shape_identity() +
    theme(legend.position= "none",
          plot.title= element_text(face="bold", size= 20, hjust= 0.5),
          axis.title= element_text(face="bold", colour="black", size= 16),
          axis.text= element_text(face="bold",size=14),
          panel.grid.minor= element_blank(),
          axis.ticks.length= unit(0.15,"cm")) +
    scale_x_continuous(expand=expansion(add= 0.1),limits=c(0, 4),
                       breaks=seq(0, 4, by=1))+
    scale_y_continuous(expand=expansion(mult= c(0.01,0.01)),limits=c(0, 120),
                       breaks=seq(0, 120, by=30))+
    scale_colour_manual(values= kovesi.rainbow(n1))
  
  p2 <-ggplot(behavior2, aes(x= Response_Time, y= Accuracy, group= Subject, 
                             color= Subject)) +  
    geom_point(aes(colour= Subject)) +
    geom_line(aes(y= my.rmc2$model$fitted.values), linetype= 1) +
    ggtitle("Nonmusic Group") +
    ylab("Accuracy (%)") +
    xlab("Response Time (s)") +
    theme_grey() +
    scale_shape_identity() +
    theme(legend.position= "none",
          plot.title= element_text(face="bold", size= 20, hjust= 0.5),
          axis.title= element_text(face="bold", colour="black", size= 16),
          axis.text= element_text(face="bold",size=14),
          panel.grid.minor= element_blank(),
          axis.ticks.length= unit(0.15,"cm")) +
    scale_x_continuous(expand=expansion(add= 0.1),limits=c(0, 4),
                       breaks=seq(0, 4, by=1))+
    scale_y_continuous(expand=expansion(mult= c(0.01,0.01)),limits=c(0, 120),
                       breaks=seq(0, 120, by=30))+
    scale_colour_manual(values = kovesi.rainbow(n2))
  
  # Combine the plots
  plot_grid(p, NULL, p1, NULL, p2, nrow= 1, rel_widths=c(1, 0.05, 1, 0.05, 1))
  
  dev.off()
  
--------------------------------------------------------------------------------  
# Calculation of balanced integration score (BIS) 
  data <- read_excel("C:/Users/ALIENWARE/Desktop/data.xlsx")
  BIS <- function(data) {
    n <- length(data$Group)    # sample size to correct var()-function result (which uses n-1)
    srt <- sqrt( ((n-1)/n) * var(data$mean_rt_c) )   # sample standard deviation across all rts
    spc <- sqrt( ((n-1)/n) * var(data$pc) )          # sample standard deviation across all rts
    mrt <- mean(data$mean_rt_c)                      # mean across all rts
    mpc <- mean(data$pc)                             # mean across all pcs
    zrt <- (data$mean_rt_c-mrt)/srt                  # standardized rts
    zpc <- (data$pc-mpc)/spc                         # z-standardized pcs
    data$bis <- zpc - zrt                            # Balanced Integration Score
    
    return(data)                                     # return data.frame with added variable 'bis'
  }
  
  data1 <- BIS(data)
 
  library(haven)
  write_sav(data1, path = "data1.sav")
    
  # install.packages("readr")
  # library(readr)
  # write_csv(data1, path = "data1.csv") 
  
--------------------------------------------------------------------------------
## Nonparametric Type III two-way mixed ANOVA on BIS using package "ARTool"
  behavior <- read_excel("C:/Users/ALIENWARE/Desktop/data-2SD.xlsx")
  names <- c('Subject','Group','Emotion')
  behavior[,names] <- lapply(behavior[,names] , factor)
  str(behavior)
  summary(behavior)
    
  m <- art(BIS~Emotion*Group+(1|Subject), data=behavior)
  summary(m)
  anova(m, response="aligned")
  b <- anova(m, response="art",type=3)
  print(b, verbose=TRUE,digits=2) 
   
  # Effect Size for two-way mixed ANOVA using package "effectsize"
  print(F_to_omega2(0.14,1,70, ci=0.95,alternative="two.sided"),digits=3)  # Group
  print(F_to_omega2(34.62,4,280,ci=0.95,alternative="two.sided"),digits=3) # Emotion
  print(F_to_omega2( 2.07,4,280,ci=0.95,alternative="two.sided"),digits=3) # Emotion*Group
  
  interpret_omega_squared(c(.000, .321, .015), rules="field2013") 
  
  # BF for main effects and interactions using "anovaBFcalc" MWS method with alpha = -0.5 (default)
  bf_mws(F=  0.14,df1=1,df2=70, report.as="BF10",alpha=-0.5) # Group
  bf_mws(F= 34.62,df1=4,df2=280,report.as="BF10",alpha=-0.5) # Emotion 
  bf_mws(F=  2.07,df1=4,df2=280,report.as="BF10",alpha=-0.5) # Emotion*Group 
  
  # Post comparison between emotion
  c <- art.con(m, ~ Emotion, adjust="bonferroni") %>%
    summary() %>% 
    mutate(sig = symnum(p.value, corr=FALSE, na=FALSE,
                        cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
                        symbols = c("***", "**", "*", ".", " ")))
  print(c, verbose=TRUE,digits=3)
  
  # Effect Size for pairwise comparison between emotion
  t_to_d(t= -0.76,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Happiness-Sadness 
  t_to_d(t=  6.38,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Happiness-Fear
  t_to_d(t=  1.06,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Happiness-Anger
  t_to_d(t= -5.22,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Happiness-Neutrality  
  t_to_d(t=  7.14,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Sadness-Fear
  t_to_d(t=  1.82,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Sadness-Anger
  t_to_d(t= -4.46,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Sadness-Neutrality
  t_to_d(t= -5.33,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Fear-Anger
  t_to_d(t=-11.60,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Fear-Neutrality
  t_to_d(t= -6.27,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Anger-Neutrality
  
  interpret_cohens_d(c(-0.05,0.38,0.06,-0.31,0.43,0.11,-0.27,-0.32,-0.69,-0.37),rules="cohen1988")
  
  # BF for paired t-test using "BayesFactor" 
  exp(ttest.tstat(t= -0.76,n1=36,rscale=0.707)[['bf']]) # Happiness-Sadness  
  exp(ttest.tstat(t=  6.38,n1=36,rscale=0.707)[['bf']]) # Happiness-Fear 
  exp(ttest.tstat(t=  1.06,n1=36,rscale=0.707)[['bf']]) # Happiness-Anger
  exp(ttest.tstat(t= -5.22,n1=36,rscale=0.707)[['bf']]) # Happiness-Neutrality  
  exp(ttest.tstat(t=  7.14,n1=36,rscale=0.707)[['bf']]) # Sadness-Fear   
  exp(ttest.tstat(t=  1.82,n1=36,rscale=0.707)[['bf']]) # Sadness-Anger  
  exp(ttest.tstat(t= -4.46,n1=36,rscale=0.707)[['bf']]) # Sadness-Neutrality
  exp(ttest.tstat(t= -5.33,n1=36,rscale=0.707)[['bf']]) # Fear-Anger  
  exp(ttest.tstat(t=-11.60,n1=36,rscale=0.707)[['bf']]) # Fear-Neutrality 
  exp(ttest.tstat(t= -6.27,n1=36,rscale=0.707)[['bf']]) # Anger-NNeutrality 
 
--------------------------------------------------------------------------------  
## Nonparametric Type III two-way mixed ANOVA on Intensity Ratings using package "ARTool"
  behavior <- read_excel("C:/Users/ALIENWARE/Desktop/data-2SD.xlsx")
  names <- c('Subject','Group','Emotion')
  behavior[,names] <- lapply(behavior[,names] , factor)
  str(behavior)
  summary(behavior)
    
  m <- art(Intensity_Ratings~Emotion*Group+(1|Subject), data=behavior)
  summary(m)
  anova(m, response="aligned")
  b <- anova(m, response="art",type=3)
  print(b, verbose=TRUE,digits=3) 
  
  # Effect Size for two-way mixed ANOVA using package "effectsize"
  print(F_to_omega2(8.47,1,70, ci=0.95,alternative="two.sided"),digits=3)  # Group
  print(F_to_omega2(74.55,4,280,ci=0.95,alternative="two.sided"),digits=3) # Emotion
  print(F_to_omega2( 6.15,4,280,ci=0.95,alternative="two.sided"),digits=3) # Emotion*Group
  
  interpret_omega_squared(c(.094, .508, .067), rules="field2013") 
  
  # BF for main effects and interactions using "anovaBFcalc" MWS method with alpha = -0.5 (default)
  bf_mws(F=  8.47,df1=1,df2=70, report.as="BF10",alpha=-0.5) # Group
  bf_mws(F= 74.55,df1=4,df2=280,report.as="BF10",alpha=-0.5) # Emotion 
  bf_mws(F=  6.15,df1=4,df2=280,report.as="BF10",alpha=-0.5) # Emotion*Group 
  
  # Post comparison between emotion
  c <- art.con(m, ~ Emotion, adjust="bonferroni") %>%
    summary() %>% 
    mutate(sig = symnum(p.value, corr=FALSE, na=FALSE,
                        cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
                        symbols = c("***", "**", "*", ".", " ")))
  print(c, verbose=TRUE,digits=2)
  
  # Effect Size for pairwise comparison between emotion
  t_to_d(t=  0.49,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Happiness-Sadness 
  t_to_d(t=  2.60,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Happiness-Fear
  t_to_d(t= -5.56,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Happiness-Anger
  t_to_d(t= 11.25,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Happiness-Neutrality  
  t_to_d(t=  2.11,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Sadness-Fear
  t_to_d(t= -6.05,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Sadness-Anger
  t_to_d(t= 10.77,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Sadness-Neutrality
  t_to_d(t= -8.16,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Fear-Anger
  t_to_d(t=  8.65,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Fear-Neutrality
  t_to_d(t= 16.81,df_error=280,paired=TRUE,ci=0.95,alternative="two.sided") # Anger-NNeutrality
  
  interpret_cohens_d(c(0.03,0.16,-0.33,0.67,0.13,-0.36,0.64,-0.49,0.52,1.00),rules="cohen1988")
  
  # BF for paired t-test using "BayesFactor" 
  exp(ttest.tstat(t=  0.49,n1=36,rscale=0.707)[['bf']]) # Happiness-Sadness  
  exp(ttest.tstat(t=  2.60,n1=36,rscale=0.707)[['bf']]) # Happiness-Fear 
  exp(ttest.tstat(t= -5.56,n1=36,rscale=0.707)[['bf']]) # Happiness-Anger
  exp(ttest.tstat(t= 11.25,n1=36,rscale=0.707)[['bf']]) # Happiness-Neutrality  
  exp(ttest.tstat(t=  2.11,n1=36,rscale=0.707)[['bf']]) # Sadness-Fear   
  exp(ttest.tstat(t= -6.05,n1=36,rscale=0.707)[['bf']]) # Sadness-Anger  
  exp(ttest.tstat(t= 10.77,n1=36,rscale=0.707)[['bf']]) # Sadness-Neutrality
  exp(ttest.tstat(t= -8.16,n1=36,rscale=0.707)[['bf']]) # Fear-Anger  
  exp(ttest.tstat(t=  8.65,n1=36,rscale=0.707)[['bf']]) # Fear-Neutrality 
  exp(ttest.tstat(t= 16.81,n1=36,rscale=0.707)[['bf']]) # Anger-NNeutrality 
  
  # Post comparison for Emotion*Group 
  c <- art.con(m, ~ Emotion*Group, adjust="bonferroni") %>%
    summary() %>% 
    mutate(sig = symnum(p.value, corr=FALSE, na=FALSE,
                        cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
                        symbols = c("***", "**", "*", ".", " ")))
  print(c, verbose=TRUE,digits=5)
  
  # Effect Size for Post comparison between music-nonmusic
  t_to_d(t= 3.51,df_error=189.57,paired=FALSE,ci=0.95,alternative="two.sided") # Happiness 
  t_to_d(t= 2.37,df_error=189.57,paired=FALSE,ci=0.95,alternative="two.sided") # Sadness 
  t_to_d(t= 2.57,df_error=189.57,paired=FALSE,ci=0.95,alternative="two.sided") # Fear 
  t_to_d(t= 2.95,df_error=189.57,paired=FALSE,ci=0.95,alternative="two.sided") # Anger 
  t_to_d(t=-0.58,df_error=189.57,paired=FALSE,ci=0.95,alternative="two.sided") # Neutrality 
  
  interpret_cohens_d(c(0.51,0.34,0.37,0.43,-0.08),rules="cohen1988")
  
  # BF for two-sample t-test using "BayesFactor" 
  exp(ttest.tstat(t= 3.51,n1=36,n2=36,rscale=0.707)[['bf']]) # Happiness 
  exp(ttest.tstat(t= 2.37,n1=36,n2=36,rscale=0.707)[['bf']]) # Sadness 
  exp(ttest.tstat(t= 2.57,n1=36,n2=36,rscale=0.707)[['bf']]) # Fear 
  exp(ttest.tstat(t= 2.95,n1=36,n2=36,rscale=0.707)[['bf']]) # Anger  
  exp(ttest.tstat(t=-0.58,n1=36,n2=36,rscale=0.707)[['bf']]) # Neutrality 
  
--------------------------------------------------------------------------------
                                  Figures
--------------------------------------------------------------------------------    
# after install the packages below, run and load them                                                                                                                                                    ,125))) 
  library(ggplot2)
  library(PupillometryR)
  library(ggpp)
  library(cowplot)
  
  # Load data  
  dat1 <- read_excel("C:/Users/ALIENWARE/Desktop/data-2SD.xlsx")
  
  # Factor variables
  dat1$Subject <- factor(dat1$Subject)
  dat1$Group <- factor(dat1$Group,levels=c("Music", "Nonmusic"))
  dat1$Emotion <- factor(dat1$Emotion,levels=c("Happiness","Sadness","Fear",
                                               "Anger","Neutrality"))
  summary(dat1)
  str(dat1)
  
## Draw Raincloud plots for Figure 2 
  tiff(filename="Fig2.tiff", width=16, height=16, units="in", res=800)
  
  # Create plot for Accuracy
  p1 <- ggplot(dat1, aes(x=Group, y=Accuracy, fill=Group, colour=Group))+
    geom_flat_violin(position=position_nudge(x=0, y=0), adjust=2, trim=FALSE)+
    geom_point(position=position_jitternudge(width=0.15, height=1, seed=123,
               x= -0.25, direction="as.is", nudge.from="jittered"), size=1.3)+
    geom_boxplot(aes(x=as.numeric(Group)+0, y=Accuracy), outlier.shape=NA,
                 alpha=0.3, width=.1, linewidth=0.5, colour="BLACK")+
    theme_test()+
    guides(fill="none", colour="none")+
    facet_wrap(~Emotion, nrow=5)+
    scale_colour_manual(values= c("#FFA60F", "#1380A1"))+
    scale_fill_manual(values= c("#FFA60F", "#1380A1"))+
    theme(axis.title=element_text(face="bold", colour="black", size=16),
          strip.text=element_text(face="bold", colour="black", size=14),
          strip.background=element_rect(color="black", fill="grey85", linetype="solid"),
          strip.placement="outside", 
          strip.switch.pad.grid=unit(7, "pt"),
          axis.text=element_text(face="plain", colour="black", size=14),
          axis.ticks=element_line(linewidth=0.3, color="black"),
          axis.ticks.length=unit(0.15,"cm"),
          panel.spacing.x=unit(1,"lines"),
          plot.title=element_text(hjust= -0.076, face="bold", size=22))+
    coord_flip()+
    labs(x="Group", y="Accuracy (%)")+
    scale_y_continuous(expand=expansion(mult= c(0.01,0.01)),
                       limits=c(-35, 140),breaks=seq(-35, 140, by=35))+
    ggtitle("A") + theme(plot.title=element_text(hjust= -0.2, vjust= 1))
  
  # Create plot for Response Time
  p2 <- ggplot(dat1, aes(x=Group, y=Response_Time, fill=Group, colour=Group))+
    geom_flat_violin(position=position_nudge(x=0, y=0), adjust=2, trim=FALSE)+
    geom_point(position=position_jitternudge(width=0.15, height=1, seed=123,
               x= -0.25, direction="as.is", nudge.from="jittered"), size=1.3)+
    geom_boxplot(aes(x= as.numeric(Group)+0, y=Response_Time), outlier.shape=NA,
                 alpha=0.3, width=.1, linewidth=0.5, colour="BLACK")+
    theme_test()+
    guides(fill="none", colour="none")+
    facet_wrap(~Emotion, nrow=5)+
    scale_colour_manual(values= c("#FFA60F", "#1380A1"))+
    scale_fill_manual(values= c("#FFA60F", "#1380A1"))+
    theme(axis.title=element_text(face="bold", colour="black", size=16),
          strip.text=element_text(face="bold", colour="black", size=14),
          strip.background=element_rect(color="black", fill="grey85", linetype="solid"),
          strip.placement="outside",
          strip.switch.pad.grid=unit(7, "pt"),
          axis.text=element_text(face="plain", colour="black", size=14),
          axis.ticks=element_line(linewidth=0.3, color="black"),
          axis.ticks.length=unit(0.15, "cm"),
          panel.spacing.x=unit(1, "lines"),
          plot.title=element_text(hjust= -0.076, face="bold", size=22))+
    coord_flip()+
    labs(x="Group", y="Response Time (s)")+
    scale_y_continuous(expand=expansion(mult= c(0.01,0.01)),
                       limits=c(-1, 5),breaks=seq(-1, 5, by=1))+
    ggtitle("B") + theme(plot.title=element_text(hjust= -0.2, vjust= 1))
  
  # Combine the plots
  plot_grid(p1, NULL, p2, nrow=1, rel_widths=c(1, 0.1, 1))
  
  dev.off() 
  
--------------------------------------------------------------------------------   
## Draw Raincloud plots for Figure 4
  tiff(filename="Fig4.tiff", width=16, height=16, units="in", res=800)
  
  # Create plot for BIS    
  p1 <- ggplot(dat1, aes(x=Group, y=BIS, fill=Group, colour=Group))+
    geom_flat_violin(position=position_nudge(x=0, y=0), adjust=2, trim=FALSE)+
    geom_point(position=position_jitternudge(width=0.15, height=1, seed=123,
               x= -0.25, direction="as.is", nudge.from="jittered"), size=1.3)+
    geom_boxplot(aes(x=as.numeric(Group)+0, y=BIS), outlier.shape=NA,
                 alpha=0.3, width=.1, linewidth=0.5, colour="BLACK")+
    theme_test()+
    guides(fill="none", colour="none")+
    facet_wrap(~Emotion, nrow=5)+
    scale_colour_manual(values= c("#FFA60F", "#1380A1"))+
    scale_fill_manual(values= c("#FFA60F", "#1380A1"))+
    theme(axis.title=element_text(face="bold", colour="black", size=16),
          strip.text=element_text(face="bold", colour="black", size=14),
          strip.background=element_rect(color="black", fill="grey85", linetype="solid"),
          strip.placement="outside", 
          strip.switch.pad.grid=unit(7, "pt"),
          axis.text=element_text(face="plain", colour="black", size=14),
          axis.ticks=element_line(linewidth=0.3, color="black"),
          axis.ticks.length=unit(0.15,"cm"),
          panel.spacing.x=unit(1,"lines"),
          plot.title=element_text(hjust= -0.076, face="bold", size=22))+
    coord_flip()+
    labs(x="Group", y="Balanced Integration Score (z)")+
    scale_y_continuous(expand=expansion(mult= c(0.01,0.01)),
                       limits=c(-9, 5),breaks=seq(-9, 5, by=2))+
    ggtitle("A") + theme(plot.title=element_text(hjust= -0.2, vjust= 1))
   
  # Create plot for Intensity Ratings
  p2 <- ggplot(dat1, aes(x=Group, y=Intensity_Ratings, fill=Group, colour=Group))+
    geom_flat_violin(position=position_nudge(x=0, y=0), adjust=2, trim=FALSE)+
    geom_point(position=position_jitternudge(width=0.15, height=1, seed=123,
               x= -0.25, direction="as.is", nudge.from="jittered"), size=1.3)+
    geom_boxplot(aes(x= as.numeric(Group)+0, y=Intensity_Ratings), outlier.shape=NA,
                 alpha=0.3, width=.1, linewidth=0.5, colour="BLACK")+
    theme_test()+
    guides(fill="none", colour="none")+
    facet_wrap(~Emotion, nrow=5)+
    scale_colour_manual(values= c("#FFA60F", "#1380A1"))+
    scale_fill_manual(values= c("#FFA60F", "#1380A1"))+
    theme(axis.title=element_text(face="bold", colour="black", size=16),
          strip.text=element_text(face="bold", colour="black", size=14),
          strip.background=element_rect(color="black", fill="grey85", linetype="solid"),
          strip.placement="outside",
          strip.switch.pad.grid=unit(7, "pt"),
          axis.text=element_text(face="plain", colour="black", size=14),
          axis.ticks=element_line(linewidth=0.3, color="black"),
          axis.ticks.length=unit(0.15, "cm"),
          panel.spacing.x=unit(1, "lines"),
          plot.title=element_text(hjust= -0.076, face="bold", size=22))+
    coord_flip()+
    labs(x="Group", y="Intensity Rating")+
    scale_y_continuous(expand=expansion(mult= c(0.01,0.01)),
                       limits=c(0, 8),breaks=seq(0, 8, by=1))+
    ggtitle("B") + theme(plot.title=element_text(hjust= -0.2, vjust= 1))
   
   # Combine the plots
   plot_grid(p1, NULL, p2, nrow=1, rel_widths=c(1, 0.1, 1))
   
   dev.off()   
     
--------------------------------------------------------------------------------